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Abstract. We show under general conditions that the linearized force-based quasi- 
continuum (QCF) operator has a positive spectrum, which is identical to the spectrum 
of the quasinonlocal quasicontinuum (QNL) operator in the case of second-neighbour 
interactions. Moreover, we establish a bound on the condition number of a matrix of 
eigenvectors that is uniform in the number of atoms and the size of the atomistic re- 
gion. These results establish the validity of and improve upon recent conjectures ([7, 
Conjecture 2] and [6, Conjecture 8]) which were based on numerical experiments. 

As immediate consequences of our results we obtain rigorous estimates for convergence 
rates of (preconditioned) GMRES algorithms, as well as a new stability estimate for the 
QCF method. 



1. Introduction 

Quasicontinuum methods are a prototypical class of multiscale models that directly 
couple multiple modeling regions to reduce the computational complexity of modelling 
large atomistic systems. These methods are useful for computing the interaction of local- 
ized material defects such as crack tips or dislocations with long-range elastic fields of a 
crystalline material. The force-based quasicontinuum (QCF) method [3, 4, 21] partitions 
the material into two disjoint regions, the atomistic region and the continuum region. 
It assigns forces to the degrees of freedom within each region using only the respective 
model, be it atomistic or continuum. This simplifies the formulation of the method as no 
special interaction rules are needed near the atomistic-continuum interface. The simplic- 
ity of mixing forces combined with the lack of spurious interface forces (so-called "ghost 
forces") make the force-based method a popular approach, and this technique is widely 
applied in the multiscale literature [1, 2, 3, 11, 13, 21, 22]. 

A potential drawback of the QCF method is that it does not derive from an energy 
(as it generally produces a non-conservative field). While the practical implications of 
this fact are still under investigation, it is already clear that the analysis of the QCF 
method poses formidable challenges. A series of recent articles has been devoted to its 
study [6, 7, 8]. For example, it was shown in [7, 8] that the linearized QCF operator is 
not positive definite, and that it is not uniformly stable (in the number of atoms and the 
size of the atomistic region) in discrete variants of most Sobolev spaces. 
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However, numerical experiments in [6, 7] showed some unexpected spectral properties. 
Conjecture 2 in [7] states that the spectrum of ^-eigenvalues of the linearized QCF opera- 
tor is identical to that of the operator associated with the quasinonlocal QC method [23] . 
This is particularly surprising since the quasinonlocal QC method is energy-based, and 
thus indicates that the QCF operator is diagonalizable and that its spectrum is real. Con- 
jecture 8 in [6] states that the condition number of a matrix of eigenvectors of the QCF 
operator grows at most logarithmically. This is an important fact for understanding the 
solution of the QCF equilibrium equations by the GMRES method. 

The purpose of the present paper is to provide rigorous proofs for these numerical 
observations. We define the QCF method and introduce the necessary notation in Section 
2. In Section 3 we establish all results for the case of next-nearest neighbour interactions 
as in the numerical experiments in [6, 7]. Then, in Section 4 we extend the results to the 
more technical case of finite-range interactions. In the case of finite-range interactions 
we cannot make the comparison between the spectra of the QCF and QNL operators. 
Instead, we prove that the spectrum of the QCF operator lies between the spectrum of 
the atomistic operator and the spectrum of the continuum operator (Theorem 12). We 
note, moreover, that we were able to construct a matrix of eigenvectors for the QCF 
operator whose condition number is bounded uniformly in the number of atoms and the 
size of the atomistic region. This result is in fact stronger than the conjectures made in 
[6]. Finally, in Section 5, we analyze variants of preconditioned QCF operators to obtain 
rigorous convergence rates for preconditioned GMRES methods as well as a new stability 
estimate. 



For the sake of brevity, we will keep the introduction to the atomistic model and the 
various flavours of quasicontinuum approximations to an absolute minimum. We refer 
to [4, 6, 7, 8, 9, 14, 15, 16] for detailed discussions. Note, in particular, that we have left 
out the usual rescaling factor e. This reduces the complexity of the notation and is justified 
since in this paper we are primarily concerned with algebraic aspects of quasicontinuum 
operators. 

2.1. Notation for difference operators. In this section, we summarize the notation 
and certain elementary results for some standard finite difference operators with periodic 
boundary conditions. 

2.1.1. Periodic domains. We identify M. N with periodic infinite sequences as follows: 



2. Formulation of the QCF Method 



R N = {u e R z : u e+N = u g for all i e Z}. 
The £ 2 -inner product on R N , and its associated norm, are defined as 



N 




We will frequently use a subspace U C Mr of mean zero functions, 



U = {ueR N : (u,e) = 0}, 
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where e = (l) £eZ e R N . 



The orthogonal projection onto U is denoted Pu : M. — >■ 

iv 



1 * 

(P w u)* = u e - — ^2 



N 

k=l 



or, m matrix notation, 



P u = I-±e®e. (1) 

2.1.2. The backward difference operator. The difference operator D : M. N — > ~R N is defined 
by 

Du t = {Du)i =ug- 

We note that rg(D) = U and ker(P) = span{e} where rg denotes the range and ker 
denotes the kernel of an operator. We also remark that, here and throughout, unless 
specifically stated otherwise, we will not distinguish between an operator and its associated 
matrix representation in M, NxN . 

2.1.3. The discrete Laplace operator. The second generic operator that we will encounter 
is the negative Laplace operator L : 1R — > WL N , 

Lu t = (Lu)e = —ui-x + 2ui - u e+1 . 

As for the difference operator, rg(L) = U and ker(L) = span{e}. 
Using summation by parts, we obtain 

(Lu, v) = (Du, Dv), 

which implies that L = D T D and hence L = L T . Since Le = 0, we have the identities 

LP U = P U L = L. (2) 

We also note that ||L|| < 4, and that this bound is attained for even N, as well as in the 
limit N — > oo. 

Since L is singular, we also define the modified negative Laplace operator 

L 1 =L + e®e = L + (I-P u ), (3) 
so that L\u = Lu if u e IA and L\e = e. This operator is invertible and satisfies 

L 1 1 L = LL 1 1 = Pu. 

2.1 A. The translation operator. The translation operator T : M. N — > ~R N is defined by 

Tu e = (Tu)e = u i+1 . (4) 

T is an orthogonal operator, i.e., T T T = I and its eigenbasis can be written explicitly as 

Tw k = X k Wk, Afc = e n , (w k ) & = e n for 1 < k < N. (5) 

The eigenvalues of T are located on the unit circle T := {i G C : \t\ — 1} and, in the 
limit — y oo, are dense in T. 

We remark that we can write the difference operator D and the negative Laplace op- 
erator as Laurent polynomials in T: D — Pd{T) and L = where 

p D (t) = (l-r 1 ), and p L (t) = (-t + 2-r 1 ). 
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In general, if p(t) is a polynomial, then the spectrum of the operator p(T) is {p(Xk) '■ 1 < 
k < N}, and the eigenvectors are the same as for T. Since T is a normal operator, all 
polynomials p(T) are also normal. 

Finally, we note that Te = e, which implies that T or any polynomial of T commutes 
with e ® e. In particular, this implies that all polynomials in T (e.g., L, D) and the 
operators L\ and Py commute. 

2.2. The linearized atomistic operators. We consider an atomistic model problem 
with periodic boundary conditions. We let U be the set of admissible displacements of 
an iV-periodic chain: the set of all iV-periodic displacements with mean zero. The latter 
condition is necessary to ensure that the systems of equations that we consider are well 
posed. If F > is a fixed macroscopic strain, then the energy (per period) of the atomistic 
chain subject to a displacement u G M. N is given by 

R N 
r=l 1=1 

where <fi G C 2 (0,+oo) is a pair interaction potential, for example, a Lennard-Jones or 
Morse potential, and R G N, R > 2, can be thought of as a discrete cutoff radius. (Note 
that, even though we have defined £ a for all u G M. N , only u G W are admitted in the 
solution of the minimization problem.) 

The Cauchy-Born or local quasicontinuum (QCL) approximation of £ a is the functional 

N N 

e=i £=i 

where W is the Cauchy-Born stored energy function, W(s) = Yl^i't'i' 1 " 3 )- 

Our analysis in the present paper concerns properties of the Hessians L a = D 2 S a (0) and 
L c = D 2 £ c (0) and the quasicontinuum operators that we derive from them. For future 
reference we write out L a and L c explicitly, 

(L & u) e = 22 Kf{- u £+t + 2ue ~ U£-r), and (6) 

r=l 
R 

(L c u) e = € F r 2 (-U£ + i + 2n £ - u^) = Wp (Lu) e , (7) 

r=l 

where the constants <^)" rF and Wp are given by 

R 

<f>';p = <j>"(rF) and W F = W"(F) = ^ <j>%r 2 . 

r=l 

We understand both L a and L c as linear operators from to Mr, defined by the above 
formulas, but are primarily interested in their properties on U. For example, we note that 
if (ftp > and Wp > 0, then both are positive definite on U and in particular invertible (see 
[7, Eq. (2.2) and Sec. 2.2] and [5, Prop. 1 and Prop. 2] for the next-nearest neighbour 
case, and [10] for finite range), however, both operators have a non-trivial kernel that 
contains e. For the continuum operator, L c , the stability condition Wp > is sharp, that 



THE SPECTRUM OF THE QCF OPERATOR 



5 



is, L c is positive definite if and only if W F > 0. We work to prove our spectral results on 
L qci up to this sharp stability criterion. 

2.3. The force-based quasicontinuum method. The force-based quasicontinuum (QCF) 
method is obtained by mixing the forces from the atomistic and the continuum model. 
To this end, we define atomistic and continuum regions A and C that satisfy 

AUC = {1,...,N} and A n C = 0. (8) 

We define the QCF forces 



_dSHv) 
qcf/ \ _ I du e ' 



Fr (u 



8£ c (u) 



if i e A, 
ifieC. 



Linearization of the nonlinear QCF operator F(u) = {F^ ci (u))f =1 at u — 0, yields the 
linear QCF operator (or simply, QCF operator), L qcf : M. N — > ~R N , 

which is the focus of our studies in the present paper. 

Unfortunately, L qcf as defined above, does not map U to U, hence we will normally 
consider the projected QCF operator (see [7, Sec. 2.3] for more detail) 

Lf = PuL*. 

To conclude this section we represent L qcf in a more compact way. By considering the 
characteristic function of A, 

_ j l, iiieA, 
Xl ~ \ 0, ttieC, 

and the associated diagonal operator X : WL N —> ~R N , 

(Xu)e = xeue, 

we can write L qcf in the form 

L qcf = [1 - X]L C + XL a = L c + X[L a - L c ]. (10) 

2.4. The quasinonlocal quasicontinuum method. The second atomistic/continuum 
hybrid scheme that will feature prominently in our investigations is the quasinonlocal 
quasicontinuum ( QNL ) method [23] . We note that the QNL method is only defined for 
second- neighbour interaction range (i.e., R = 2). Extensions to further neighbours ex- 
ist [9, 12, 20], but we only use the version up to second neighbours in this paper. The 
QNL method is conservative with energy functional 



eeAuc eeA 



eec 
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The linearized QNL operator is the Hessian of £ qnl at u = 0, that is, L qnl = D 2 S qnl (0). 
The operator L qnl : M. N — > M^, is most easily written in variational form [5, Sec. 3.3], 

(L qnl u, u) = ^ Du e Dvi - $ F Lu e Lv e . (11) 
eeAuc eeA 

Based on this representation one can show that, if W F > and <p F > 0, then L qnl is 
positive definite on U (see [5, Prop. 3] for the case when $jf < 0; the case 4>2F > is 
trivial) . 

3. The £ 2 -Spectrum of the Second-Neighbour L qcf Operator 

In [7] the invertibility of the QCF operator was investigated analytically and numer- 
ically, and several interesting numerical observations were left as conjectures. Similar 
observations were also used in [6] to study the performance of iterative solution methods 
for the QCF operator. In the present section we rigorously establish some of the most 
fundamental of these conjectures in the next-nearest neighbour case. We will then extend 
the results, to finite range interactions in Section 4. 

3.1. Similarity of L qcf and L qnl . In [6, 7] it was observed in numerical experiments 
that the spectra of L qcf and L qnl coincide. In this section we provide a rigorous proof by 
explicitly constructing a similarity transformation between L qcf and L qnl . The main ideas, 
after which the proof will be straightforward, are the following two new representations 
of the L qcf and L qnl operators. 

Proposition 1. Let R = 2, then L qcf and L qnl have, respectively, the representations 

L qcf = W' F L - <p2 F XL 2 , and (12) 
L qnl = WpL - <Pl F LXL. (13) 



Proof. We begin by noting that, for R = 2, the operators L a and L° may be written as 

L a = <P F L + ^ F [AL - L 2 ] = W'^L - <f>% F L 2 , 
L c = 4>' F L + ^ F [4L] = W£L. 

Using these formulas, the operator L qcf (as defined in (9)) can be written in terms of 
the atomistic and the continuum operators (10) as 

L q C f = XL a + [/ _ x]L c = X [W'^L - <pl F L 2 ] +[I-X] [W'^L] . 

From this we immediately obtain (12). 

To rewrite the QNL operator, we note that we can write (11) as 

(L qnl u, v) = W F {Du, Dv) - ^ F {XLu, Lv) 
= W F (D T Du, v) - $ F (LXLu, v), 

for all u, v G M. N , and we therefore obtain (13). □ 
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Based on (12) and (13) we will deduce the similarity of the QCF and QNL operators. 
Since L is not invertible, we introduce the nonsingular operator L\ : M. N — > W N defined by 
(3). The following result confirms Conjecture 2 in [7] (a related conjecture with different 
boundary conditions is Conjecture 6 in [6]). 

Theorem 2 (Similarity of L qcf and L qnl ). If R = 2 then the operators Lf and L qnl 
are similar, with similarity transformation L\ defined by (3): 

Lf = L^U^U. 

In particular, the spectra of Lf and L qnl coincide. 

Proof. Using formulas (2), (12), and (13), a straightforward computation yields the desired 
identity: 

L x Lf = L 1 P u L qci = LL qci 
= W F LL - ^pLXL 2 

= [WpL - (j)2 F LXL\Li = L qnl L l . □ 

3.2. Condition number of the £ 2 eigenbasis. Since L qnl is self-adjoint, there exists 
an orthonormal matrix V qnl G M> NxN , and a diagonal matrix A containing the eigenvalues 
of L qnl , such that 

L qn V qnl = V qnl A. (14) 

Note in particular that A also contains the zero eigenvalue. Since the operators L qnl and 
L qcf are similar, there exists also an invertible operator V qcf G ]R ArxAr such that 

L qcl V qci = V qci A. 

As suggested by Theorem 2, a possible choice for the eigenvectors is L^~ 1 l /qnl , since in 
that case we have 

Lf L^ l V qnl = L^ l L qnl V qnl = LfV qnl A. 

However, cond(L^ 1 l /qnl ) = cond(Lj~ 1 ) = 0(N 2 ), which is much worse than the numerical 
observations in [7] and suggests a poor scaling of the eigenvectors. 

To produce a better eigenbasis, it is important to note that the choice of eigenvectors 
is not unique even after fixing the ordering as we are always free to rescale them. This 
turns out to be a crucial ingredient in our following construction of an eigenbasis with a 
uniformly bounded condition number. The following result is inspired by Figure 4.2 in 
[7] and Conjecture 8 in [6]. It does not precisely confirm these, but establishes a closely 
related and in fact stronger result. 

Theorem 3. Suppose that R = 2, then the operator V qcf : M. N — > M. N , 

V qci = [W';i - ^ F P u XL]V qn \ (15) 

diagonalizes Lf , that is, Lfv qcf = V qci A, where A is the diagonal matrix of eigenvalues 
associated with L qnl (14). Moreover, if Wp > and cf)" F > then cond(U qcf ) is bounded 
above by a constant that depends on 4>2p/Wp, but is independent of N and A. 
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Remark 1 . The choice of V qcf is motivated by the following calculation. Starting with 
the similarity result of Theorem 2, we derive 

Lf = L^L^Li = L^V qnl A(V qnl ) T L 1 , 

and we scale the eigenvectors by A, which gives 

L -lyqnl A = ^-lyqnl^yqnljTyqnl = ^-l^qnly qnl 

= L^[W F L - ^ F LXL}V^ 1 = [W'^Pu - <& F PuXL]V**. 

This is equivalent to the choice of V qci in (15) when restricted to U. □ 

Proof of Theorem 3. Step 1: Diagonalization. In a straightforward computation we 
obtain 

Lfv^= [w;l-^ f p u xl 2 ] [w^i-^ f p u xl]v^ 

= [Wp - (% F P U XL\ [W F L - <P'i F LXL] V qnl 
= [W F I - (j)'i F PuXL\ V^ l A 
= V^A. 

Step 2: Estimating cond(V qcf ). We now assume that W F > and cj)" F > 0. To estimate 
cond(V qcf ) we can ignore the positive constant multiple W F as well as the orthonormal 
matrix V qnl , that is, we have 

cond(y qcf ) = cond(A) = ||A|| HA" 1 ]], (16) 
where A = I — aPyXL, 

with constant a = t^r, and the convention H^l" 1 !! = +oo if A is not invertible. We note 

that the condition <ft F , W F > implies that a < 1/4. 
Elementary estimates give the upper bound 

\\A\\ < 1 + H||P W ||||X||||L|| = 1 +A\a\. (17) 

We similarly get the lower bound 

\\Au\\ > 1- |a|||P w ||||X||||L|| = 1 — 4|afj for all = 1, (18) 

which gives an estimate for ||^4 _1 || whenever W F + 402 F > 0. In the following we prove 
a bound for ||^4 _1 || that holds whenever W F > 0, that is, up to the sharp stability limit, 
which is a more involved result. 

To estimate we use the facts that (i) HA" 1 1| = HA" 2 " || , and (ii) if 

\\A T u\\ > 7o||u|| VuGl^, 

for some constant 7 > 0, then A T is invertible and ||^4~ r || < l/7o- In Lemma 4 we 
establish precisely this fact, assuming that a < 1/4, with a constant 70 that depends only 
on a but not on N or A. □ 

A generalization of the following technical lemma used in the previous proof will also be 
required in the finite range interaction case. It follows from Lemma 11 by choosing Z = I. 



THE SPECTRUM OF THE QCF OPERATOR 9 

Lemma 4. Let a < 1/4, then there exists a constant 70 > 0, which depends on a but is 
independent of N and of A, such that 



I - aP u XL} T u\\ > 7o||u|| Vu E 



Remark 2. From the proof of Lemma 11 we see that, in the case (fr'^p < 0, the constant 
7o is explicitly given by 

7 2 = 1 + 8a 2 - Wa 2 + 4a 4 , (19) 
where a = (p F /W F , and the resulting condition number estimate by 

cond(V qcf ;W) 2 < (1+4| a|)2 =: c(a) 2 . 

1 + 8a 2 - AVa 2 + 4a 4 

If 4>2 F is moderate but W F — > 0, then a — > —00. A brief calculation shows that in this 
limit c(a) behaves asymptotically like 

c(a) ~ 2 5/ V + 0(|a| 3/2 ) as a ->■ 00. □ 



4. Finite Range Interactions 

In this section we extend the results of Section 3 to the case of finite range interactions 
(i.e., with arbitrary finite R). We begin by stating a simplified formulation of our main 
results. We note, however, that our actual results are more general than the following 
theorem. In particular, we can replace the assumption §" rF < 0, r > 2, by a more general 
condition. This theorem will be proved in Section 4.4. 

Theorem 5. Suppose that <f>" rF < for 2 < r < R and that cf)" RF 7^ 0. 

(i) There exists an operator V qcf : M. N — > M, N , which diagonalizes Lq C{ , that is, 

LfV qci = V qci A, 

where A is a diagonal real matrix of eigenvalues (Xj)^ =l . 

(ii) If Wp > then V q is invertible and cond(l/ qcf ) is bounded above by a constant 
that depends on the coefficients 4>" F , r = 1, . . . , R, but is independent of N and A. 

(Hi) IfWp > and if the eigenvalues are ordered, then 

A • < A, < A*, 

where (X C j)^ =1 and (X^)j =1 denote the ordered eigenvalues of, respectively, L c and L a . In 
particular, we have the bounds 

X 1 = 0, and AW F sin 2 (f ) < Xj < A<j)" F for j = 2, . . . , N. 
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4.1. Symmetrization of L qcf . We recall from (10) the definition of the finite range QCF 
operators, 

L q C f = l° + X[L*-L c ], and Lf = L c + P u X[L a - L% 

where L a and L c are, respectively, the atomistic and the continuum operators defined in 
(6) and (7). In terms of the translation operator defined in (4), we can express L a and L c 
as 

R R 

L c = J2r 2 <j) / ; F [-T + 2I -T- 1 ], and L a = £ & F [ - T r + 21 - T~ r ] . (20) 

r=l r=l 

To show that LQ Cf is diagonalizable, we aim to construct a matrix Y such that YL qci Y _1 
is symmetric. We first note that we can write the difference L a — L c as a Laurent poly- 
nomial of the translation operator T: 

R 

L a -L c = €f [(-T r + 21- T~ r ) - r 2 (-T + 21- T' 1 )] 

r=2 

= E €f [(T r - I) (T~ r -I)-r\T- I) (T^ 1 - /)] . 

r=2 

Thus, if we define the Laurent polynomial 

R 

Kt) = £ € F [(t r - m~ r - 1) - r\t - l)(r J - 1)] for t G C \ {0}, (21) 

r=2 

then we obtain L a — L c = 6(T). The crucial observation is the following: If we can 
factorize b(t) as b(t) = p(t)p(l/t) then we can choose Y = p(T) to obtain 

L a -L c = b(T) = piT^T" 1 ) = p(T)p(T T ) = YY T = Y T Y, (22) 

and we immediately obtain 

YL qci = YL C + YX [Y T Y] = [L c + YXY T ] Y. (23) 

Here, we have used the fact that all operators that are polynomials in T commute. A 
similar result holds also for L[} cf , though this would require further preparation. If Y were 
invertible (this will turn out to be false), then this would imply that L qcf is similar to a 
symmetric, hence, normal matrix, and is therefore diagonalizable. 

The desired polynomial factorization result is essentially the Riesz-Fejer factorization 
lemma [18, Sec. 53], which we state and prove in a slightly more general form. 

Lemma 6. Let b(t) be a Laurent polynomial with real coefficients such that b(t) = b(l/t). 
Then there exists a polynomial pif) such that bit) = p(t)p(l/t). 

If, in addition, b(t) > for all t G T := {t G C : \t\ = 1}, then p(t) can be chosen to 
have real coefficients. 



Remark 3. In Lemma 6, if b(t) < on T, then it can be factorized as b(t) = —p(t)p(l/t), 
where p(t) has real coefficients. □ 
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Proof of Lemma 6. Let at~ R be the leading term in b(t) with a^O, then a(t) := t R b(t) 
is a polynomial with a(0) 7^ 0. Moreover, a and b share the same roots, which we collect 
into a set A, so that 

b(t) = ar R Y[(t-\) m{x \ 

AeA 

where m(A) denotes the multiplicity of A. 
Next, we define the auxiliary polynomial 

Pl ( t )= n (t-xr^. 

AeA,|A|<l 

Since b(t) has real coefficients it follows that A G A if and only if A G A, and m(A) = m(A), 
and hence pi(t) also has real coefficients. We also note that 

pi(i/o= n (-A/tr (A) (t-i/A) m(A) . 

AeA,|A|<l 

Since b(t) = b(l/t), for each root A with |A| > 1, 1/A is also a root, with the same 
multiplicity, and therefore we have 

b(t)=p 1 (t)p 1 (l/t) h{t), 

where 

6 1 (t) = « 1 r fcl (t-A) m(A) 

AeA,|A|=l 

for some constants «i 6 K and k\ G N. By construction, b\ has real coefficients and 
= 6i(l/*)- 

We will later assume that b has no roots on T except for t = 1, and for that case, the 
proof would be complete. To establish the lemma in its full generality, we now distinguish 
two cases. 

Case 1: b(t) > on T ■ For t G T, we have 1/t = i, and hence, upon reordering the 
product in 



P i(t) P1 (i/t) = n [(t - x) m ^\t- A) m w 



AGA,|A|<1 



> for all t G T. 



Thus, if > on T, then we also have b\ (t) > on T . This implies that the multiplicity 
m(A) of all roots A G T fl A is even. We can therefore define 

P2(t)= n ^-Ar (A)/2 , 

AeA,|A|=l 

and argue similarly as above, to obtain that p2{t) has real coefficients and that 

bi(t) = a 2 r k *p 2 (t)p 2 (l/t) =: b 2 (t)p 2 (t)p 2 (l/t), 

for some constants a 2 G M and k 2 G N, and b 2 (t) = a 2 t~ h2 . We can again deduce that 
b 2 (t) = b 2 (l/t) and that b 2 (t) > on T, which implies that k 2 = and a 2 > 0. Thus, we 
obtain p{t) by defining 

P(t) = y/aipi(t)p2(t). 
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Case 2: b(t) ^ on T ■ In this case the multiplicity of roots A G T may be odd, and 
therefore we define 

P2(t)= n (t-xr^. 

AeA,|A|=l,ImA>0 

Thus, P2{t) contains all roots of b\(t) with positive imaginary part and p 2 (l/t) contains 
all roots of b\(t) with negative imaginary part. We are only left to find any roots at ±1. 
Since bi(t) = bi(l/t) it follows that these roots must be even (possibly m(— 1) = 0), and 
hence we can define 

p 3 (*) = (t- l)m(l)/2( t+1 )m(-l)/ 2j 

to obtain 

h(t) = a 2 t k *p 2 (t)p 2 (l/t)p 3 (t)p 3 (l/t) 

for some constants a 2 G M. and k 2 G Z. Arguing as before we find that a 2 > and k 2 = 0, 
and hence we obtain the result if we define 

Pit) = V^Pl(t)P2(t)p 3 (t). □ 

For any Laurent polynomial b(t) we call a polynomial p(t) that satisfies b(t) = p(t)p(l/t) 
a GRF-factorization (generalized Riesz-Fejer factorization) of b(t). 

Since the Laurent polynomial b(t) defined in (21) clearly satisfies the condition b(t) = 
b(l/t), we have now indeed established (22) and (23). Thus, if Y were invertible, then we 
could deduce that L qcf (and similarly Lg Cf ) is diagonalizable. Unfortunately, as follows 
from the next lemma, Y is always singular. 

Lemma 7. Let p(t) be any GRF-factorization of the Laurent polynomial b(t) defined in 
(21), then there exists a polynomial pi(t) such that pit) = it — l) 2 p\(t). In particular, the 
operator Y = p(T) can be rewritten as 

Y = LY X where Y x = -T Pl (T). 



Proof. It can be immediately seen from the definition of b(t) that b(l) = and, since 
b(l) = p(l) 2 , it follows that p(l) = 0. Next, dividing b(t) by (t - l)^ -1 - 1) we obtain 

P(t)p(l/t) b(t) 

(f-l)(t-i-l) 

R 

= [(! + t + • • • + < r_1 )(l + ^ + ' • • + 

r=2 

and hence 

p(t)p(l/t) 



r+lN _ r 21 



(f-l)(f-l-l) 



= 0. 



Thus, we see that 1 is also a (multiple) root of ^z^fprzji ' anc ^ we can conc hide that 
p(t) = (t - 1) 2 PiW for some polynomial pi (t) . 

Upon rewriting the representation p(t) = (t — l) 2 p\(t) as 

P {t) = [(t - lXr 1 - 1)] [ - tp x (t)] = [ - 1 + 2 - r 1 } [ - t Pl (t)] 

we obtain Y = p(t) = LY\. □ 
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Thus, as a consequence of Lemma 7, we see that Y is always singular. Nevertheless, in 
order to construct a similarity transformation to diagonalize -^o° f , we can make a similar 
modification as in the next-nearest neighbour case, and simply replace L by its invertible 
variant L x defined in (3). Of course this still leaves us to verify that Y x is invertible, for 
which we will introduce conditions in Sections 4.2 and 4.4. 



Proposition 8. Let p(t) be any GRF-factorization ofb(t), and let Y = p(T) = LY\, as 
in Lemma 7, then 

[LiY 1 ]I%*=[L° + YXY T \[L 1 Y 1 ]. 
In particular, ifY\ is invertible then Lq C{ is diagonalizable. 

Proof. We use (22) to represent L a — L c , the fact that polynomials of T and L\ commute, 
and the fact that PuL\ = L x Pu — L, to obtain 

[L x Y x ]Lf = [L X Y X ]L C + Y X [L X P U }XLY X Y?L 
= L C [L X Y X ] + Y X LXLY X T P U [L X Y X } 

= [L C + Y 1 LXLY X T ][L 1 Y X ]. □ 

4.2. Condition number of the £ 2 eigenbasis. To investigate the invertibility of Y x we 
will study the Laurent polynomial 

bl M = b M (24) 

[(t-iX^-i)] 2 

on the unit circle T. It will quickly become apparent that b x {t) needs to be bounded 
away from zero on T in order to obtain invertibility of Y x and bounds on the inverse 
that are uniform in N and A. Consequently, we focus on the case where b(t) does not 
change sign on T (note that b(t) is always real on T). To show that this is a reasonable 
assumption we will, in Section 4.4, study the case of non-positive coefficients cj)" F < for 
r = 2, . . . , R, which is of particular interest in applications as this condition is satisfied 
for most practical interaction potentials. We will show that one can obtain bounds of the 
following type: there exist positive constants f3o, f3 x such that 

^o<\bi(t)\<p 2 x VtGT. (25) 

Clearly, since b(t) is real on T, it is necessary that b x (t) does not change sign on T, that 
is, either b x (t) > or b x (t) < 0. If b x (t) > then the matrix Y x defined in Lemma 7 is 
real, however, if b x (t) < then it is imaginary. For the following analysis we prefer to 
construct a real coordinate transform. 

According to Remark 3 and Lemma 7, we can choose a polynomial p(t) = (t — l) 2 p x (t) 
with real coefficients such that 

b(t) = ap(t)p(l/t), 

where a G {+1, —1} is the sign of b x (t) on T. Thus, if we (re-)define Y x accordingly as 

Y x = -Tp x (T) (26) 
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then Proposition 8 implies that 

[L^Lf = [L c + aiLY^XiLY^] [L&]. (27) 

Next, we show that (24) implies invertibility of Y"i, including explicit bounds. Beforehand, 
however, we make a brief remark on the connection of the sign of b(t) and a relationship 
between L a and L c . 

Remark 4- The assumption b\(t) > on T, or, more generally, b(t) > on T, implies 
that 

(L a u,u)>(L c u,u) WueR N . 

This follows simply from the fact that L a and L c share the same eigenvectors, and that 
the spectrum of L a — L c , which is contained in b(T) is non-negative. 
Conversely, if b(t) < on T then, by the same argument, 

(L a u,u) < (L c u,u) VhgM" 

These observations are interesting in their own right, and will moreover lead to useful 
estimates on the eigenvalues of the QCF operator L^ 1 in Section 4.3. □ 

Lemma 9. Let Y\ be defined as in (26) and suppose that b(t) satisfies (25), then Y\ is 
invertible and we have the bounds 

IN|</3!, and || Yf 1 !! < l//3„. 



Proof. Let p(t) — (t — l) 2 pi(t) be a GRF-factorization of o~b(t) with real coefficients, then 
Pi{t) is a GRF-factorization of with real coefficients. Hence, for t G T, we obtain 

b(t)= Pl (t) Pl (t) = \ Pl (t)\ 2 , 

which implies that 

A)<|tpi(t)|<ft VtGT. (28) 
By assumption, Y\ = — Tpi(T) =: qi(T). Since T is unitary, this implies that Y\ 
has an orthonormal eigenbasis and that the spectrum of Y\ is contained in qi(T). This 
immediately gives 

|| YJ < SU P |tpi(t)| < 0!. 

teT 

Moreover, we have 

HFfl-^ inf ||y lU || >inf \t Pl (t)\ >p , 

||w||=l tST 

which gives the second bound. □ 

According to (27), and Lemma 9, we have 

[LyY^Lf [LxY]" 1 = L c + aiLYjXiLY^ =: L s ^, (29) 

that is, L qcf and L sym are similar matrices. Since L sym is real and symmetric, there exists 
an orthogonal operator V^ sym G ]R ArxAr such that 
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where A is the diagonal matrix of eigenvalues of L-Q Cf , and we obtain 

L qci = \yf 1 Li 1 V sym ]A\YC 1 L^ 1 V Bym ]- 1 , 

that is, the operator Y i ~ 1 Li 1 V aym diagonalizes L qci . As in the nearest neighbour case, the 
eigenbasis Y^~ L x V sym is poorly scaled and would lead to an 0(N 2 ) condition number. 
However, the same argument as in Remark 1 shows how to rescale the basis to obtain the 
following theorem. 

Theorem 10. Suppose that the Laurent polynomial b{t) defined in (21) satisfies (25), 
and let Y\ be defined by (26). Then the operator V qc{ G M ArxAr , 

V qc{ = [WpYf 1 + aP u XY^L] V sym 

diagonalizes L^ ci , that is, Lo cf V qcf = ^ qcf A, where A is a real diagonal matrix of eigen- 
values. 

Moreover, if Wp > 0, and if a/3f/Wp > —1/4, then V qcf is invertible and cond(l/ qcf ) 
is bounded above by a constant that depends only on Wp,(3o,(3i, and, in particular is 
independent of N and A. 

Proof. The proof of this result is very similar to the proof of Theorem 3, and hence we 
shall be fairly brief. First, we note that, by Lemma 9, the matrix Y\ can be chosen to be 
invertible, and hence V qci is well-defined. Moreover, we recall that Y\ is a polynomial in T 
and hence commutes with all operators that commute with T, such as other polynomials 
in T, the modified negative Laplace operator L 1; and the projection operator Py. 

Step 1: Diagonalization. As in the computation in the proof of Theorem 3, we obtain 

L qc V qcf = [W F L + aPuXY^L 2 ] [W'jYf 1 + aP u XY?L] V sym 

= [WpY^ 1 + aP u XY?L] [W£L + oLY x P u XYfL\ V sym 

= [WpYf 1 + aP u XY^L\ V sym A 

= V qci A. 

Step 2: Estimating cond(y qcf ). Suppose now that Wp > 0. As before, estimating 
||\^ qcf || is straightforward. Using Lemma 9, we obtain 

\\V qcf || < W^IIYT 1 !! + \\P U XY^L\\ < W'p[/f3 + 4ft, (30) 

where (5q and 0i are the constants from (25). 

To estimate V q from below, we will use again the fact that || |V qcf ] _1 || = ||[U qcf ]~ T || < 
l/(Wp%) where % > satisfies 

||[y qcf ] T «|| >To|H| Vn G M. N . (31) 
Writing out (V qci ) T in full, we get 

[V qcf f = W£Yr T [V 8ym f[l + ^Y^LY.XPu] 
= W' F [V* ym } T Y- T [I - aY^LYiXPu] , 
where a = —ojW'p. If a/3 2 < 1/4, then we can use Lemma 11 below to deduce that 
\\[l+^nY^LY l XP u \v\\ >7o|MI V^l^, 
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where 70 depends only on a/3 2 = —a(3f/Wp, but is independent of N and A. This implies 
(31) with 70 = W^ 7o /ft. 

Combining (30) and (31) gives the stated result. □ 

Lemma 11. Let A = I — aPuXZ T LZ , where Z G IR ArxAr commutes with Pu, and where 
a£l satisfies 

-00 < a\\Z\\ 2 < 1/4, 
then there exists a constant 70 > 0, depending only on a\\Z\\ 2 , such that 

\\A T u\\ > 7o||w|| Vjig1 n 



Proof. We decompose A T into 

A T = I + aZ T LZXP u — [I — P u ] +Pu[l~ aZ T LZXP u ] , 

where we have used the fact that P^L = L and that Pu commutes with Z. Since Pu is 
an orthogonal projection we obtain, again using PuL = L, 

\\A T v\\ 2 = - P u }v\\ 2 + \\[I - aZ T LZX}P u v\\ 2 . (32) 

We will show next that 

|| [1 -aZ T LZX]w\\ 2 > (1 -e)\\w\\ 2 Vu> G U, (33) 

where e G (0, 1) is defined in (37) and depends only on a||Z|| 2 , but not on N or A. Hence, 
(33) combined with (32), gives the desired result 

\\A T v\\ 2 > min(l, 1 - e)(||[l - P u ]v\\ 2 + ||P^;|| 2 ) = 7o ||i;|| 2 , 



with 7 = a/1 — e- 

Proof of (33). We begin by splitting the operator, using X 2 = X, into 
[I - aZ T LZX] = X[I- aZ T LZX] + [I - X] [I - aZ T LZX] 
= X[I - aZ T LZ]X + [I -X][I - aZ T LZX] 
=: Si + S 2 . 

Since X is an orthogonal projection, we have, for any wsW, 

\\[I -aZ T LZX}w\\ 2 = \\Snu\l 2 +\\S 2 w\\ 2 . (34) 

Estimating S\. Since Si = X[I — aZ T LZ]X is a symmetric operator, the following 
variational bound is found using the fact that < (Lx, x) < 4. We obtain 

(X[I - aZ T LZ]Xw, w) = \\Xw\\ 2 - a(LZXw, ZXw) 

> min (1,1 -4a||Z|| 2 )||Xw||, 

where the last equality follows by distinguishing the cases a < and a > 0. 
In summary, if ct||Z|| 2 < 1/4, then we have the ^^.-independent bound 

||Siw|| 2 > [min(l,l -4a||Z|| 2 )] 2 ||Xw|| 2 Vw G R N . (35) 
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Estimating S 2 . Due to the good estimate on Si we only need fairly rough estimates on 
the term HS^u* || 2 . Application of the Cauchy-Schwartz inequality and a weighted Cauchy 
inequality provides the estimate 

\\S 2 w\\ 2 > (1 - e)\\[I - X]w\\ 2 + (1 - e-VHtJ - X][Z T LZ]Xw\\ 2 , 

for any e G (0, 1). Using the fact that / — X is an orthogonal projection, ||L|| < 4, we can 
further estimate 

\\S 2 w\\ 2 > (1 - e)||[J - X]w\\ 2 + (1 - e-^iellZllVUX^II 2 . (36) 

Combining the estimates. If we define a = Aa\\Z\\ 2 , and insert (35) and (36) into (34), 
we obtain, for all w G 14, 

|| [/ - aZ T LZX]w\\ 2 > { min(l, (1 - a) 2 ) + (1 - e" 1 )^} ||X«;|| 2 + (1 - e) || [/ - X]w\\ 2 

> min { min(l, (1 - a) 2 ) + (1 - e _1 )d 2 , 1 - e} |H| 2 , 

for any e G (0, 1). Since min(l, (1 — a) 2 ) > it is clear that choosing e sufficiently close 
to 1 gives a positive lower bound. To optimize this constant with respect to e, we have to 
choose e to equalize the two competitors in the outer min formula. The resulting choice 
is 

d 2 + ±d 4 -|d 2 , d<0, 

+ v/2(a 2 - a 3 ) + a 4 , < d < 1, 

which concludes the proof of (33). (As a matter of interest, e — > 1 as a — > 1, e = for 
a — 0, and e — >• 1 as 5 — >• —00.) □ 

4.3. Estimates for the eigenvalues. Using the similarity to a symmetric matrix that 
we have established in the previous section we can now give sharp bounds on the spectrum 
of Lf . 

Theorem 12. Suppose that (25) holds, then LQ Ct has a real, ordered spectrum (Xj)^ =1 . 
If we denote the ordered eigenvalues of L a and L c , respectively, by (A*)^ =1 and (X C j)^ =1 
then 

Aj c < Xj < AJ, for j = 1, ... , N, if b(t) > 0, and 
X^XjKX*, forj = l,...,N, ifb(t)<0. 



a — a 2 



Proof. We know from Theorem 10 that Lq c1 is diagonalizable and that it is similar to the 
self-adjoint operator L sym defined in (29), which has a real spectrum that is identical to 
the spectrum of Lq C{ . We will next show that, for all u G M. N , 

(L c u, u, ) < (L sym w, u) < (L a w, u), if a = 1, and 

(L a w, u, ) < (L sym M , u) < (L c u, u), if a = -1. ^ 

From these inequalities, the min- max characterisation of eigenvalues [17, Sec. XIII. 1] 
immediately gives the stated result. 

To prove (38) we will take the following starting point: 

(L sym u, u) = (L c u, u) + a(Y T XYu, u) = (L c u, u) + a(XY T u, XY T u). (39) 
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From here on, we treat the cases o = 1 and o = — 1 separately. 

Case 1: a = 1. Using (39) and the fact that Y T Y = L a — L c , we have 

(L sym u, u) < (L c u, u) + (Y T u, Y T u) = (L c u, u) + ([L a - L c ]u, u) = (L a u, u). 

For the lower bound we use (39) and the fact that (XY T u,XY T u) is non-negative to 
obtain 

(L sym u,u) > (L c u,u). 

Case 2: a = — 1. The idea for a = — 1 is essentially that the roles of L a and L c are 
reversed. Note that, now, Y T Y = L c — L a . Hence, for the upper bound, we get 

(L sym u,u) < (L c u,u), 

whereas, for the lower bound, we get 

(L sym u, u) > (L c u, u) - (Y T u, Y T u) = (L c u, u) - ([L c - L a ]n, u) = (L a n, u). □ 

4.4. The case of non-positive coefficients. For many of the common interaction po- 
tentials, such as the Lennard- Jones potential, <f>(r) = Ar~ 12 + Br~ 6 , or the Morse poten- 
tial, 4>{r) = exp(— 2a(r — r )) + 2exp(— a(r — r )), there exists a minimal strain F* such 
that, for all F > F*, we have 

4>" F < for r = 2, . . . , R. (40) 

In most cases, it is reasonable to assume that the macroscopic strain F remains in this 
region, as it would require extreme compressive forces to violate it. Hence, for the re- 
mainder of the section, we will assume that (40) is satisfied. The following two lemmas 
reduce this case to the one studied earlier in this section. 

Lemma 13. Suppose that (40) holds, then b(t) > in T. 

Proof. Similarly as in the proof of Lemma 7 we rewrite b(t) in the form 

m = X>:v[(* r - m~ r - 1) - r\t - lx*- 1 - 1)] 

r=2 

R 

= (t- i)(r 1 €f [(* r_1 + t r ~ 2 + ■■■ + i)(^ r+1 + t~ r+2 + . . . + 1) - r 2 ] . 

r=2 

(41) 

It is easy to see that (t — l)^" 1 — 1) is non-negative on T, and moreover, for t £ 7", 

(f 1 + f~ 2 + ... + i)(r r+1 + r r+2 + . . . + 1) = If" 1 + t r " 2 + . . . + 1| 2 < r 2 . (42) 
Hence we obtain the stated result. □ 

Lemma 14. Suppose that (40) holds and that <fi" RF < 0. Then (25) is satisfied with 
constants /3q and fi\ that are independent of N and A. 
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Proof. The upper bound 0i in (40) obviously exists since bi (t) is a continuous function 
on a compact set T. 

To show the existence of the lower bound it is sufficient to show that 

&i(t)>0 for all teT. (43) 

Suppose that < at some point t\ G T. Then b(t\) < and hence at least one 

term in the definition oib(t) (21) is non-positive. To be precise, there exists r G {2, . . . , i?} 
such that <j)" F < and 

i-€ F ) [r\h - i)(tr x - 1) - (*; - mr - 1)] < o 

It follows from (41) and (42) that this may only happen at t\ — 1. However, a striaght- 
forward computation shows that 

= lim ^-i)(r'-i)-y-i)(r--i) = ^ 

Hence no point t\ G T such that &i(ti) < exists. □ 
We are now in a position to complete the proof of Theorem 5. 

Proof of Theorem 5. Item (i) is a special case of Proposition 8, taking into account that, 
for non-positive coefficients, b(t) defined in (21) is non-negative and Y\ can therefore be 
chosen to be real. Item (ii) follows from Theorem 10 and Proposition 14. Item (iii) is 
established in Theorem 12. □ 

We conclude this section with a result that gives the sharp bounds on /3 , fi\ for the case 
of non-positive coefficients. 

Proposition 15. Let b\(t) be defined by (24) and suppose that (40) holds; then (25) 
holds with constants 

= Y,{-€ F ) 2r2 + { ~ iy ~\ and 



„ r 2 (r 2 — 1) 



12 

The lower bound is attained at t = —1 and the upper bound is attained at t = 1. 
The proof of this proposition is based on the following technical lemma. 

Lemma 16. For r G N, r > 2, let f r : (0, |] ->■ R be defined as 

frW := (r 2 - ^) , (44) 

sm p \ sm p J 
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then 



inf f r ((3)=f r (n/2) = r 2 - 1 ( ^ and (45) 

0<|j3|<§ ^ 

sup / r (/3)=lim/ r (/3) = ir 2 (r 2 -l). (46) 

0<|/8|<§ 



Proof. Proof of (45). First, notice that 



/ r (7r/2) = r — sin — = r 



m _ 2 l-(-l) r 



2 



which proves that the left-hand side of (45) is not larger than the right-hand side. To 
prove the that f r (/3) > / r (7r/2) for all /3, transform 



2 sin rp \ 1 / 2 2 • 2 a sin r ^ 

1 2 = 2 r ~ r S111 P 2 

sin f3 V sin P J sin P \ sin /3 

1 / 9 9 „ sin 2 r/3\ 9 1 /r 2 sin 2 2/3 sin 2 r/3\ 

Vcos 2 /3- — +r 2 = -^- — 2 f - — + r 2 (47) 



sin 2 /3 sin 2 j3 ) sin 2 /3 \ 4 sin 2 /3 sin 2 (3 

r 2 / sin 2 2/3 sin 2 r(3 



sin /3 \ 4 r 



2 



+ r 2 



and consider the three cases: 0</3<^, ^</3<|-^;, and f — ^ < P < 1. 
Case i. (0 < /3 < I 1 ) Further transform the function / r (/3) in (47): 

/r(/3) = ( sinc2 ^ - sinc2 ^ + r '- 

sin p 

The expression in the brackets is positive since sincx = is a decreasing function for 
< x < 7r/2. This proves f r (/3) > r 2 > / r (?r/2). 

Case 2. (f- < /3 < - - ^) In this case sin 2(3 > hence 



2r — " — r 2r ' " — r ' 

1 /r 2 sin 2 2/3 sin 2 r(3\ 

fr(P) = —2—o A ■ —a ^TIT + r 

sm p \ 4 sm [3 sm p 



1/1 sin 2 r/3\ 
> —2^ - — TTT +r 2 > r 2 > / r (7r/2). 



sin p \sin /3 sin (3 

Case 3. (| — £ < /3 < 1) Denote x = | — /? (0 < x < 27) and consider the two different 
subcases: r is even/odd. 

Case 3.1. (r is even) In this case sin (3 = cosx, sin 2/3 = sin2x, and sin 2 r/3 = sin 2 (rx). 
Hence f r (9) as expressed in (47) takes the form 

„._. r 2 /sin 2 2x sin 2 rx\ 9 

/>•(/?) = — —7 -H+ r 



cos 4 x \ 4 r 2 

(sine 2 2x — sine 2 rx) + r 2 > r 2 = f r {jx j2\ 



r 2 x 2 



cos 4 x 
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Case 3.2. (r is odd) In this case sin (3 = cosx, sin 2 2(3 = sin 2 2x, and sin 2 r(3 = cos 2 (rx). 
Hence (47) transforms into 



COS 4 X 



sin 2x cos 2 rx 



> r J 



1 1 

2 + 2 



4 r 2 

2 



r 2 - 1 



+ r 2 > r 2 



/ r (7T/2). 



COS X 



cosrx 



1 + cos 2x 
2 cos rx 



Here we used the fact that 1 > cos2x > cosrx (0 < x < f-)- 
Proof of (46). First compute 

efos v r 2 sin 2 /3 - sin 2 r/3 , r 2 (/3 2 - /3 4 /3) - (r 2 (3 2 - r 4 /3 4 /3) + 0(/3 6 ) 
hm / r (p) = hm — — = hm 



P-Xi 



/3->0 



sin {3 



p->o 



— + — = -r 2 (r 2 - 1). 



To prove the inequality f r (/3) < |r 2 (r 2 — 1) consider the two cases: ^ < (3 < | and 
< /? < 7. 



Case 1. (- < (3 < \ ) In this case (46) follows directly from the following computation: 



l r \r 2 -l)^)]' 1 >\r\r 2 -V) 



sin (3 



7T 



-1 



-(r 2 - 1) sin 2 (3 
3 



> - r z - 1 sin^ - > - r z - 1 - 



3 V ' r ~ 3 
Case 2. (0 < (3 < We need to prove 
1 



4r 2 - 4 4 1 



3r 2 



> = 1. 

3 3 



MP) - ^ 2 (r 2 - 1) 



^ - ^(r 2 - 1) < (V/3 G (0, 7r/r)) (48) 



for integer r, but instead we prove that it is valid for all real values of r G [2, 00). 
First, notice that the following calculation 



fM 



- r 

3 V 



r=2 



1 1 

4 sin 2 f3 



4 



sin 2 2/3 
sin 2 (3 



4 sin 2 /3 



(4-4 cos 2 0) - 1 = 



proves (48) for r = 2. As is shown below, 

^ {fM 1 



<9r 



r 2 -!' 



< 



(V/3e(0,7r/r)), 



(49) 



which concludes the proof of (48). 

Proof of (49) . We use the standard inequalities 



sinc'(x) > — x/3, Vx > 
cosx < sinc 3 (x) Vx G [0, 7r] 
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and obtain 



and hence 



d ff r (P)\ d ( 1 sin 2 r/3\ d //3 2 sinc 2 r/3 



dr \ r 2 J dr \sin 2 /3 r 2 sin 4 /3/ dr \ sin 4 /3 
2/3 3 sine r/3 sine' rj3 ^ 2/3 3 sine r/3 (r/3) 
sin 4 /3 3 sin 4 (3 



d ffr(8) I, 2 A 2/3 3 sine r/3 (r/3) 2r 2r //3 4 sine r/3 
r - 1 < 



9r ^ r 2 3 / 3sin 4 /3 3 3 \ sin 4 /3 

3 V sm 4 /3 / 3 \ snr/3 



Proof of Proposition 15. To reduce the problem to the statement of Lemma 16 we make 
the substitution (3 = e 2tf3 , —ir/2 < 9 < tt/2, which gives 

_ ^ xll N r 2 (e 2 ^ - l)(e~ 2 ^ - 1) - (e 2i ^ - l)(e" 2 ^ - 1) 



r=2 



[(e 2 ^ - l)(e- 2 ^ - 1 



A r 2 4 sin 2 /3- 4 sin 2 r/3 _ A / r (/3) 

Z — ' 16 sin p 4 

r=2 r=2 

Application of Lemma 16 gives Proposition 15. □ 



5. Analysis of Preconditioned Lq Operators 

In this final section we present two further interesting applications of our foregoing 
analysis. First, we discuss the GMRES solution of a linearized QCF system. We rigorously 
establish some conjectures used in [6] and briefly discuss their consequences. Second, we 
prove a new stability result for the linearized QCF operator in a discrete Sobolev norm, 
which we hope will become a useful tool for future analyses of the QCF method. 

We assume throughout this section that (28) holds and that Y\ is defined by (26). More- 
over, we recall the definition of L sym from (29). Since the results are fairly straightforward 
corollaries from our analysis in Section 4 we will derive them in a less formal manner. 

5.1. GMRES-Solution of the QCF system. We consider the linearized QCF system 

Lfu = f, (50) 

where / £ U, which is to be solved for u £ U. If this system is solved using the GMRES 
algorithm (see [19] for a general introduction and [6] for a detailed discussion of using GM- 
RES for solving the QCF system), then standard estimates on GMRES convergence [19], 
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along with the analysis of the previous sections, show that the residual of the m-th iterate, 
r O) — j _ L,^ c u^ m \ satisfies the bound 

(i I — \ m 
- V7 ) ||rW||, 
1 + ViJ 

where 7 = A2/Aat = 0(1/N 2 ) (see also [6, Prop. 22]). The fraction A2/Aat is used instead 
of X\/Xn since we are solving the system in U. This convergence is rather slow and hence 
two variants of preconditioned GMRES (P-GMRES) algorithms were suggested in [6]. 
The preconditioner used was the negative Laplacian L. We will use the preconditioner Li 
instead of L, however, this is purely for the sake of a consistent notation since L^L^ = 
L _1 LQ Cf (note that rgLQ Cf = U and that L" 1 is well-defined on U). 

The first variant of P-GMRES that was considered in [6] was the standard left-precon- 
ditioned GMRES algorithm where GMRES is applied to the preconditioned system 

L^Lfu = L?f. (51) 

To obtain convergence rates, we require bounds on the eigenvalues and eigenbasis of 
L^Lf (see [6, Sec. 6.2]). 

The second variant considers again the left-preconditioned system (51) but this time 
the residual is minimized in the norm induced by the operator L\. The convergence rates 
of the resulting method are governed by the spectrum and eigenbasis of the operator 
L~ 1/2 LfL- l/2 (see [6, Sec. 6.3]). 

5.1.1. Diagonalization. We consider L^ 1 ^ 2 L^ 1 ' 2 first. Using (29), and the fact that 

— 1/2 

L± commutes with L{Yi, we obtain 

[L 1 Y 1 ][Li 1/2 L'fLi 1/2 ][L 1 Y 1 ]- 1 = L~ 1/2 L sym L~ 1/2 = V sym A(V Sirm ) T , 

where A is the real diagonal matrix of eigenvalues and V sym an orthonormal matrix of 
eigenvectors of L^ 1 ^ 2 L sym L\ 1 ^ 2 . Hence, we conclude that L^^L^L^ 1 is diagonalizable 
with real spectrum A: 

[L~ l/2 Lf L~ l/2 ][Y^ l L^V sym ] = [Y^L^V^jA. (52) 

— 1/2 

Multiplying the equation by L x , we obtain 

[L^Lf^L-^V 3 ^] = [Yf^-^V^JA, (53) 

thus showing that also L^Lq is diagonalizable with the same real spectrum A. We note 
that this rigorously establishes a variant of [6, Conjecture 10]. 

5.1.2. Condition number bounds. Using the fact that U sym is orthogonal, and Lemma 9 
to bound cond(Fi) < Pi/Po, we can obtain the following upper bounds on the condition 
number of the matrices of eigenvectors: 

cond(yr 1 Lr 1 t> sym ) < iV 2 /V/3o, and (54) 
cond(Ff 1 L- 3/2 U sym ) < iV 3/2 /V/3 . ( 55 ) 
This rigorously establishes [6, Conjectures 12 and 13]. 
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Since these bounds are not uniform in N the question arises whether we can define a 
better scaling for the eigenvectors to improve them. Note, however, that [7, Thm. 4.3] 
implies that cond(LQ Cf ) > N 1 ^ 2 , and hence no choice of eigenbasis can achieve an upper 
bound in (54) that is uniform in N. Moreover, numerical experiments in [6, Sec. 3, Figs. 
2 and 3] indicate that our bounds may be optimal. 

5.1.3. Eigenvalue bounds. To establish convergence rates for the P-GMRES solution of the 
QCF system, we still need to obtain bounds on the eigenvalues contained in A. Let (A n )^ =1 
denote the ordered eigenvalues of A, and let (A^)^ =1 and (A^)^ =1 denote, respectively, the 
ordered eigenvalues of L a and L c . Since L c = WpL, we know that 

XI = 0, and \ c n = Wp for n = 2, . . . , N. 

— 1/2 

In view of Remark 4, replacing k by L 2 u in the formulas, we obtain that 

either <[L~ 1/2 L a L7 1/2 ]u, u) > {[L~ l/2 L c L~ 1/2 ]u, u) Vu G U, 

or <[L^ 1/2 L a L7 1/2 ]u, u) < ([L~ 1/2 L c L~ 1/2 ]u, u) Vu G U. 

Hence, we can repeat the proof of Theorem 12 verbatim to show that 

min(A a , X c n ) <X n < max(A a , X c n ) for n = 1, . . . , N. 

At this point we need to make an assumption on the stability of the atomistic system. 
We assume that the macroscopic strain F is chosen so that 

CoIKH 2 < (L a u,u) < cju'll 2 VuGW. (56) 

The upper bound can be obtained in a straightforward computation that gives a constant 
Ci depending only on the coefficients <fr" F , r = 1, . . . , R. The lower bound means that the 
homogeneous deformation Fx lies in the region of stability of the atomistic energy (see 
[5] for a detailed discussion of this point, in particular, that c is indeed independent of 
N). For example, in the case of non-positive coefficients Remark 4 shows that this bound 
holds with c = Wp, and that in the case R = 2 one can choose c = mm(Wp, <f> F ). 

Upon noting that the stability assumption (56) is equivalent to the statement that 
A2 > c and A^ < c\, we can now deduce that 

min(c , Wp) <X n < max(ci, Wp), 
which are bounds that are independent of and A. 

5.1.4. Convergence rates for P-GMRES. From the foregoing discussion we obtain the 
following convergence rates for the P-GMRES solution of (50) (see [6, Sec. 6.2 and Sec. 
6.3] for details of these derivations): 

For the standard left-preconditioned GMRES algorithm we obtain bounds on the pre- 
conditioned residual, 

||Lr> M || < CA^ 3 g m ||L7 1 r (0) ||, (57) 
where C > and q G (0, 1) are independent of A^ and A. 

For the left-preconditioned P-GMRES algorithm, which minimizes the preconditioned 
residual in the norm induced by Li, we obtain 

|| L -i/2 r ( m .)|| < CN 2 q m \\L~ 1/2 r {0) \\, (58) 
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where C > and q G (0, 1) are independent of N and A. 

We also note that a finer analysis (see [6, Sec. 6.2 and 6.3]) shows that both variants 
of P-GMRES reduce the residual to zero in at most 0(#A) iterations. 

5.2. Stability of L^ ct in IA 2,2 . We define the discrete Sobolev-type norm 

|H| 2 = H-MI for « € E w , 

which is a norm on the space U of mean-zero functions and denote the space U equipped 
with || • || 2 by U 2,2 . Moreover, we denote the space IA equipped with the norm || • || =: || • || 
by IA ' 2 . We are interested in the question whether : IA 2,2 —> IA 0,2 is stable, uniformly 
in iV and A. 

To begin with, we note that 

H/rqcfwiM-i . , Il^o^llo _ ■ f ll^o Cf ^rVllo . 

IK ° } h ^ 2 ^~uA } \\Lu\\ ~fA} H/llo ' 

thus, the question reduces to the analysis of the operator Lq C L± l . Using the representa- 
tion 

Lf = L C + aP u X(L a - L c ) = W F L + aP u XY^Y x L 2 , 
where Y\ is defined in (26), we obtain 

Lfh- X x = W'^P U + aP u XY?Y x L. 

We now argue similar as in the proof of Theorem 9. On the space U the operator WpPu + 
PuXY^YiL can be replaced by 

W£l + aPuXYfYtL. 

For this modified operator Lemma 11 shows that it is invertible and provides uniform 
bounds on the inverse. Restricting the argument back to U we obtain the following 
theorem. 



Theorem 17. Suppose that Wp > 0, that (25) holds, and that a(3f/Wp > -1/4; then 
Lg Cf is invertible and ||(LQ Cf ) _1 || i( - w o,2 jW 2,2) is bounded above by a constant that depends only 
on Wp,l3o,l3i (that is, on the coefficients ^>" rF ^r = 1, . . .,R) but is independent of N and 
of the choice of A. 

6. Conclusion 

We have established a comprehensive £ 2 -theory of a linearized force-based quasicon- 
tinuum (QCF) operator making several conjectures from previous work [6, 7] regarding 
its spectrum and eigenbasis rigorous (at least up to a modification of the boundary con- 
ditions). We have given elementary derivations in the case of next-nearest neighbour 
interactions but have also provided proofs for arbitrary finite range interactions. Finally, 
as an immediate corollary of our analysis we have also obtained a new stability estimate 
in the space IA 2 ' 2 . 

Our results heavily use the fact that the nonlinear QCF operator is linearized at a 
homogeneous deformation and a question of immediate relevance is whether our results 
can be generalized, at least partially, to linearizations around non-uniform states. Even 
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in small neighbourhoods of homogeneous deformations it is unclear whether this can be 
done. 

Secondly, a generalization to 2D or 3D would have immense consequences as no ap- 
proach to the analysis of QCF in 2D or 3D exists at this point. 
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